#include "image.h"

//constructor
image::image() {
 
  ifstream myimage;
  vector<double> temp;
  string input;
  int t;
  double s;
  epsilon = 0.0001; //arbitrary

  cout << "Please enter input file: " << endl;
  cin >> input;
  myimage.open(input.c_str());
  
  cout << "Please enter cutoff  value: " << endl;
  cin >> s;
  cutoff = s;
  //cout << "your cutoff value is: " << cutoff << endl;

  if (myimage.is_open()) {
    while(myimage.good()) {
      for (int i=0;i<2;i++) {
	myimage >> t; //gets the integers describing size
	size.push_back(t);	
      }
      temp.clear();
      for (int i=2;i<size[0]+2;i++) { //gets the doubles describing the image
	temp.clear();
	for (int j=0;j<size[1];j++) {
	  myimage >> s;
	  temp.push_back(s);
	}
	initial.push_back(temp);
      }
    }
    myimage.close();
  }
  else {
    cout << "Unable to open image." << endl; //error message
  }
  
}

//blank deconstructor
image::~image() {

}

//displays a chosen matrix of given size
void image::print(vector < vector <double> > matrix, int m, int n) {
  if (!matrix.empty()) {      //displays the first m rows and the first n columns, even if matrix is larger
    for(int i=0;i<m;i++) {
      for(int j=0;j<n;j++) {
	cout << fixed << setprecision(0) << matrix[i][j] << "   ";
      }
      cout << endl;
    }
  }
  else {
    cout << "No file has been read." << endl;
  }
}

//transposes a matrix of any given size and returns the matrix
vector <vector <double> > image::computetranspose(vector < vector <double> > matrix, int m, int n) {

  vector <double> temp;
  vector < vector <double> > tempmatrix;
  double t;
  temp.clear();
  tempmatrix.clear();
  
  for(int i=0;i<n;i++) {
    for(int j=0;j<m;j++) {
      t = matrix[j][i];
      temp.push_back(t);
    }
    tempmatrix.push_back(temp);
    temp.clear();
  }

  return tempmatrix;

}


//multiplies two matrices and returns their product
vector < vector <double> > image::multiply(vector < vector <double> > matrix1, vector < vector <double> > matrix2, int m, int n, int p, int q) { //multiplies matrix 1 (mxn) by matrix 2 (pxq)

  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  double t;
  tempvector.clear();
  tempmatrix.clear();

  if(n==p) {
    for(int i=0;i<m;i++) {
      tempvector.clear();
      for(int j=0;j<q;j++) {
	t = 0;
	for(int k=0;k<n;k++) {
	  t = t + matrix1[i][k]*matrix2[k][j];
	}
	tempvector.push_back(t);
      }     
      tempmatrix.push_back(tempvector);
    }
  }
  else { // Returns error if matrices are of incompatible size
    cout << "Cannot multiply matrices, incompatible size." << endl;
  }
  return tempmatrix; // returns a matrix of size mxq
}

//subtracts two matrices entry by entry and returns their differences
vector < vector <double> > image::subtract(vector < vector <double> > matrix1, vector < vector <double> > matrix2, int m, int n, int p, int q) { //matrix 1 minus matrix 2

  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  double t;
  tempvector.clear();
  tempmatrix.clear();

  if(m == p && n ==q) {
    for(int i=0;i<m;i++) {
      tempvector.clear();
      for(int j=0;j<n;j++) {
	t = 0;
	t = matrix1[i][j]-matrix2[i][j];
	tempvector.push_back(t);
      }
      tempmatrix.push_back(tempvector);
    }
  }
  else { // Returns error if matrices are of incompatible size
    cout << "Cannot subtract matrices, incompatible size." << endl;
  }
  return tempmatrix;
}

//checks to see if a matrix is smaller than epsilon entry-wise. Returns true if they are, returns false if not
bool image::check_tolerance(vector < vector <double> > matrix, int m, int j) { 

//checks to see if each entry is within epsilon of 0

  for(int i=0;i<m;i++) {
    if (matrix[i][j] > epsilon || matrix[i][j] < epsilon*(-1)) {
      return false;
    } 
  }
  return true;
}


//creates a random basis of doubles all less than 1 //RAC change this to make largest entry == 1
void image::createRandomBasis(int size) {

  double temp = 0;
  vector <double> tempvector;

  srand(unsigned(time(NULL)));

  for(int i=0;i<size;i++) {
    tempvector.clear();
    for(int j=0;j<size;j++) {
      temp = ((double)rand()/(double)RAND_MAX);
      tempvector.push_back(temp);
    }
    randomBasis.push_back(tempvector);
  }
}

int image::computeEigenmatrix(vector < vector <double> > v, int size) {

  //Takes in an nxn matrix and computes its eigenvectors. We assume the matrix is square, so only one size argument is necessary.

  vector < vector <double> > matrix1;
  vector < vector <double> > matrix2;
  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  double tempvalue = 0;
  int counter = 0;
  tempvector.clear();

  for(int i=0;i<size;i++) { //fill matrix1 with random basis
    for(int j=0;j<size;j++) {
      tempvalue = randomBasis[i][j];
      tempvector.push_back(tempvalue);
    }
    matrix1.push_back(tempvector);
    tempvector.clear();
  }
  
  //fill matrix 2 with v * matrix 1, then gram schmidt it
  matrix2 = multiply(v,matrix1,size,size,size,size);
  matrix2 = gramSchmidt(multiply(v,matrix1,size,size,size,size),size,size);

  for(int i=0;i<size;i++) { //Program works one column at a time to find eigenvectors.

    tempvalue = 0;
    counter = 0;
    
    do {
      
      matrix1 = matrix2;
      matrix2 = gramSchmidt(multiply(v,matrix2,size,size,size,size),size,size);
      //print(matrix2,size,size);
      //cout << endl;
      
      tempmatrix.clear();
      tempmatrix.push_back(computetranspose(matrix2,size,size)[i]);
      tempmatrix = computetranspose(tempmatrix,1,size);
      tempvalue = computeEigenvalue(v,tempmatrix,size);
      //Calls computeEigenvalue (below) to calculate the corresponding eigenvalue for the eigenvector we just found.
      counter=counter+1;
      //cout << "counter: " << counter << endl;
    } while(check_tolerance(subtract(matrix1,matrix2,size,size,size,size),size,i)==0 && counter < 50);
    
    if(i>0 && counter<50 && (tempvalue > cutoff && tempvalue < eigenvalues[i-1])) {
      //cout << "eigenvalue is : " << tempvalue << endl;
      eigenvalues.push_back(tempvalue);
      eigenvectors.push_back(computetranspose(matrix2,size,size)[i]);
      //cout << "the eigenvector is:" << endl;
      for(int k=0;k<size;k++) {
	//cout << eigenvectors[i][k] << endl;
      }
    }     
    

    if (i == 0 && counter<50 && tempvalue > cutoff) { //RAC change epsilon to some set "cut off" value?
      //cout << "eigenvalue is : " << tempvalue << endl;
      eigenvalues.push_back(tempvalue);
      eigenvectors.push_back(computetranspose(matrix2,size,size)[i]);
      //cout << "the eigenvector is:" << endl;
      for(int k=0;k<size;k++) {
	//cout << eigenvectors[i][k] << endl;
      }
    }

    //If eigenvalue is below the decided tolerance level, it is discarded.
    //Loop does not proceed, because all future eigenvalues will be below tolerance.

    if (i > 0  && (tempvalue <= cutoff || tempvalue > eigenvalues[i-1])) {
      //cout << "found bad eigenvalue" << endl;
      return i;
    }

    if (i == 0 && tempvalue <= cutoff) {
      //cout << "found bad eigenvalue" << endl;
      return i;
    }

    if (counter == 50) { // If the multiplication fails to converge, we assume that the eigenvalue is small
      //cout << "LOOPED TOO MUCH!!" << endl;
      return i;
    }

  }

  return size; //Returns the number of 'good' eigenvalues. All others are less than epsilon and can be discarded.
}

double image::computeEigenvalue(vector < vector <double> > matrix, vector < vector <double> > v, int size) {

  vector < vector <double> > tempmatrix;
  double eigenvalue;

  //cout << "the column is: " << endl;
  //print(v,size,1);
  tempmatrix = multiply(matrix,v,size,size,size,1);
  //cout << "A times v is: " << endl;
  //print(tempmatrix,size,1);

  for(int i=0;i<size;i++) {
    if(tempmatrix[i][0]!=0 && v[i][0]!=0) {
      eigenvalue = tempmatrix[i][0] / v[i][0];
      if (eigenvalue != 0 && !(eigenvalue != eigenvalue)) {
	return eigenvalue;
      }
    }
  }

  return 0;

}


//gram schmidts a matrix and return the GS'd and normalized matrix
vector < vector <double> > image::gramSchmidt(vector < vector <double> > x, int size1, int size2) { //size 1 = rows, size 2 = columns

  vector <double> tempvector1;
  vector <double> tempvector2;
  tempvector1.clear();
  tempvector2.clear();
  double tempvalue=0;
  double tempvalue1=0;
  double tempvalue2=0;

  for(int i=0;i<size2;i++) {

    tempvector1.clear();
    tempvector2.clear();

    for(int j=0;j<i;j++) {

      tempvalue1=0;
      tempvalue2=0;

      for(int k=0;k<size1;k++) {
	tempvalue1 = tempvalue1 + x[k][j]*x[k][i]; //dot product of jth column with ith column
	tempvalue2 = tempvalue2 + x[k][j]*x[k][j]; //dot product of jth column with jth column
      }

      //cout << endl;
      //cout << j << " dot product " << i << " is: " << tempvalue1 << endl;
      //cout << j << " dot product" << j << " is: " << tempvalue2 << endl;
      //cout << endl;
      tempvector1.push_back(tempvalue1); //holds dot product of i with j-1...0 columns
      tempvector2.push_back(tempvalue2); //holds dot product of j-l with j-l....0 with 0 columns
    }

    //print(x,size1,size2);

    for(int j=i-1;j>=0;j--) {
      tempvalue = tempvector1[j]/tempvector2[j]; //divides dot product of jth column with ith columns by jth column with jth column
      //cout << "ratio of dot product " << j << " with " << i << " by " << j << " with " << j << " is: " << tempvalue << endl;

      for(int k=0;k<size1;k++) {
	x[k][i] = x[k][i] - tempvalue*x[k][j];
	//cout << "x[" << k << "][" << i << "] is: " << x[k][i] << endl;
      }

      //print(x,size1,size2);

    }

  }

  //cout << endl; 
  //cout << "before normalization: " << endl;
  //print(x,size1,size2);
  
  //normalize
  
  for(int i=0;i<size2;i++) {
    tempvalue1 = 0;
    for(int k=0;k<size1;k++) {
      tempvalue1 = tempvalue1 + x[k][i]*x[k][i]; //computes i.i
    }
    //cout << i << "." << i << " is: " << tempvalue1 << endl;
    tempvalue1 = sqrt(tempvalue1); //compute ||i.i||
    for(int k=0;k<size1;k++) {
      x[k][i] = x[k][i]/tempvalue1;
    }
  }
  //cout << "the matrix is: " << endl;
  
  return x;
}

void image::computeSingularValues() {

  double temp;
  
  for(unsigned int i=0;i<eigenvalues.size();i++) {
    temp = sqrt(eigenvalues[i]);
    singularvalues.push_back(temp);
  }
  return;
}

vector < vector <double> > image::scalarMultiplyColumn(vector < vector <double> > m, double d, int column, int size1, int size2) { //column is the column # you want to multiply, size1 is the number of rows, size2 is the number of columns

  double temp = 0;
  vector <double> tempvector;
  vector < vector <double> > tempmatrix;
  tempvector.clear();
  tempmatrix.clear();

  if(column < size2) {
    for(int i=0;i<size1;i++) {
      tempvector.clear();
      for(int j=0;j<size2;j++) {
	if(j == column) {
	  temp = d * m[i][column];
	  tempvector.push_back(temp);
	}
	if(j!= column) {
	  temp = m[i][j];
	  tempvector.push_back(temp);
	}
      }
      tempmatrix.push_back(tempvector);
    }
  }
  else if (column >= size2 ){
    cout << "column choice not in range" << endl;
  }
  return tempmatrix;

}

void image::computeEigenrange() {

  vector < vector <double> > tempmatrix;
  
  //cout << "initial:" << endl;
  //print(initial,size[0],size[1]);
  //cout << "eigenvectors:" << endl;
  eigenvectors = computetranspose(eigenvectors,eigenvalues.size(),size[1]);
  //print(eigenvectors,size[1],eigenvalues.size());


  tempmatrix = multiply(initial,eigenvectors,size[0],size[1],size[1],eigenvalues.size());

  //cout << "initial times eigenvectors is: " << endl;
  //print(tempmatrix,size[0],eigenvalues.size());

  for(unsigned int i=0;i<eigenvalues.size();i++) {
    tempmatrix = scalarMultiplyColumn(tempmatrix,1./singularvalues[i],i,size[0],eigenvalues.size());
  }
  
  eigenrange = tempmatrix;

  //cout << "product matrix times singular values by column is: " << endl;
  //print(eigenrange,size[0],eigenvalues.size());

  return;
}

void image::computeSVD() {

  vector < vector <double> > singular;
  vector < vector <double> > eigentrans;
  vector < vector <double> > tempmatrix;
  double temp;
  vector <double> tempvector;
  singular.clear();
  tempvector.clear();
  tempmatrix.clear();

  for(unsigned int i=0;i<eigenvalues.size();i++) {
    tempvector.clear();
    for(unsigned int j=0;j<eigenvalues.size();j++ ) {
      temp = 0;
      if (i == j) {
	temp = 1;
      }
      if (i != j) {
	temp = 0;
      }
      tempvector.push_back(temp);
    }
    singular.push_back(tempvector);
  }
  
  for(unsigned int i=0;i<eigenvalues.size();i++) {
    singular = scalarMultiplyColumn(singular,singularvalues[i],i,eigenvalues.size(),eigenvalues.size());
  }

  eigentrans = computetranspose(eigenvectors,size[1],eigenvalues.size());

  //cout << "the W matrix is:" << endl;
  //print(eigenrange,size[0],eigenvalues.size());
  //cout << "the sigma matrix is:" << endl;
  //print(singular,eigenvalues.size(),eigenvalues.size());
  //cout << "the V matrix is:" << endl;
  //print(eigentrans,eigenvalues.size(),size[1]);
  
  tempmatrix = multiply(eigenrange,singular,size[0],eigenvalues.size(),eigenvalues.size(),eigenvalues.size());
  tempmatrix = multiply(tempmatrix,eigentrans,size[0],eigenvalues.size(),eigenvalues.size(),size[1]);
  SVD = tempmatrix;

  cout << "the SVD is:" << endl;
  print(SVD,size[0],size[1]);
  
  return;
}
